#include <xtiffio.h>
#include <geo_normalize.h>

#include "OutFile64.h"

#include <stdio.h>
#include <iostream>
#include <math.h>
#include <genutils.h>
#include <timeutils.h>
#include <nc4utils.h>
#include <string>
#include <float.h>

#include <boost/algorithm/string/trim.hpp>
#include <boost/algorithm/string/case_conv.hpp>
#include <boost/algorithm/string/predicate.hpp>

#include <geo_tiffp.h>

using namespace std;

//------------------------------------------------------------------------------
// OutFile::ProductStuff
//------------------------------------------------------------------------------

OutFile::ProductStuff::ProductStuff(int32_t width,
        const productInfo_t* productInfo) {
    this->width = width;
    this->productInfo = allocateProductInfo();
    copyProductInfo(this->productInfo, productInfo);
    dataStorage = UByteDS;
    scaleType = Linear;
    scale = 1.0;
    offset = 0.0;
    minOutputVal = 0;
    maxOutputVal = 255;
    minVal = minOutputVal;
    maxVal = maxOutputVal;
    missingValue = fillPix;
    lineData = (double*) allocateMemory(width * sizeof (double),
            "OutFile::ProductStuff::lineData");
}

OutFile::ProductStuff::ProductStuff(const OutFile::ProductStuff& pStuff) {
    width = pStuff.width;
    productInfo = allocateProductInfo();
    copyProductInfo(productInfo, pStuff.productInfo);
    dataStorage = pStuff.dataStorage;
    scaleType = pStuff.scaleType;
    scale = pStuff.scale;
    offset = pStuff.offset;
    minOutputVal = pStuff.minOutputVal;
    maxOutputVal = pStuff.maxOutputVal;
    minVal = pStuff.minVal;
    maxVal = pStuff.maxVal;
    missingValue = pStuff.missingValue;
    lineData = (double*) allocateMemory(width * sizeof (double),
            "OutFile::ProductStuff::lineData");
}

OutFile::ProductStuff::~ProductStuff() {
    freeProductInfo(productInfo);
    free(lineData);
}

/**
 * set the scale factors.  note that default minOutputVal=0, maxOutputVal=255
 * @param min min geophysical value
 * @param max max geophysical value
 * @param log do you want log10 scaling
 */
void OutFile::ProductStuff::setScale(double min, double max,
        ScaleType scaleType) {
    this->scaleType = scaleType;
    minVal = min;
    maxVal = max;
    switch (scaleType) {
    case Linear:
        offset = min;
        scale = (max - offset) / (maxOutputVal - minOutputVal);
        break;
    case Log:
        offset = log10(min);
        scale = (log10(max) - offset) / (maxOutputVal - minOutputVal);
        break;
    case ArcTan:
        offset = min;
        scale = max;
        minVal = calcPhysicalVal(minOutputVal);
        maxVal = calcPhysicalVal(maxOutputVal);
        break;
    default:
        printf("-E- OutFile::setScale - invalid scaleType = %d\n",
                (int) scaleType);
        exit(EXIT_FAILURE);
    }
}

void OutFile::ProductStuff::setScale(double min, double max,
        ScaleType scaleType, double minOutput, double maxOutput) {
    minOutputVal = minOutput;
    maxOutputVal = maxOutput;
    setScale(min, max, scaleType);
}

/**
 * set the scale factors.  note that default minOutputVal=0, maxOutputVal=255
 * @param scale slope
 * @param offset intercept
 * @param scaleType type of scaling to calculate
 */
void OutFile::ProductStuff::setScaleOffset(double scale, double offset,
        ScaleType scaleType) {
    this->scaleType = scaleType;
    this->scale = scale;
    this->offset = offset;

    // have to set these so calcPhysicalVal does not limit the physical val
    minVal = 0 - FLT_MAX;
    maxVal = FLT_MAX;
    minVal = calcPhysicalVal(minOutputVal);
    maxVal = calcPhysicalVal(maxOutputVal);
}

void OutFile::ProductStuff::setScaleOffset(double scale, double offset,
        ScaleType scaleType, double minOutput, double maxOutput) {
    minOutputVal = minOutput;
    maxOutputVal = maxOutput;
    setScaleOffset(scale, offset, scaleType);
}

double OutFile::ProductStuff::calcOutputVal(double val) const {
    double outVal;

    if (val == badPixelValue)
        return missingValue;

    if (val == landPixelValue)
        return landPix;

    // don't scale if output type is floating point
    if (dataStorage == FloatDS || dataStorage == DoubleDS)
        return val;

    switch (scaleType) {
    case Linear:
        outVal = (val - offset) / scale;
        break;
    case Log:
        if (val < 0)
            return minOutputVal;
        else
            outVal = (log10(val) - offset) / scale;
        break;
    case ArcTan:
        outVal = scale * ((atan(0.5 * val - offset) / atan(offset)) + 1);
        break;
    default:
        printf("-E- OutFile::ProductStuff::calcOutputVal - invalid scaleType = %d\n",
                (int) scaleType);
        exit(EXIT_FAILURE);
    }

    if (outVal < minOutputVal)
        return minOutputVal;
    if (outVal > maxOutputVal)
        return maxOutputVal;
    return outVal;
}

double OutFile::ProductStuff::calcPhysicalVal(double val) const {
    double physicalVal;

    if (val == missingValue)
        return badPixelValue;

    // don't scale if out output type is floating point
    if (dataStorage == FloatDS || dataStorage == DoubleDS)
        return val;

    switch (scaleType) {
    case Linear:
        physicalVal = val * scale + offset;
        break;
    case Log:
        physicalVal = pow(10, val * scale + offset);
        break;
    case ArcTan:
        physicalVal = (tan((val / scale - 1) * atan(offset)) + offset) / 0.5;
        break;
    default:
        printf("-E- OutFile::calcPhysicalVal - invalid scaleType = %d\n",
                (int) scaleType);
        exit(EXIT_FAILURE);
    }

    if (physicalVal < minVal)
        return minVal;
    if (physicalVal > maxVal)
        return maxVal;
    return physicalVal;
}

//------------------------------------------------------------------------------
// OutFile
//------------------------------------------------------------------------------

OutFile::OutFile() {
    width = 0;
    height = 0;
    qualityData = NULL;
    currentLine = 0;
    colorType = Grayscale;
    fileMinVal = DBL_MAX;
    fileMaxVal = 0 - DBL_MAX;
    resolution = 0;
    deflate = 0;

    red = (uint8_t*) allocateMemory(256, "red");
    green = (uint8_t*) allocateMemory(256, "green");
    blue = (uint8_t*) allocateMemory(256, "blue");

    // init to grayscale
    for (int i = 0; i < 256; i++) {
        red[i] = i;
        green[i] = i;
        blue[i] = i;
    }

    metaData = (meta_l3bType*) allocateMemory(sizeof (meta_l3bType),
            "OutFile::metaData");
    metaData->north = 90.0;
    metaData->south = -90.0;
    metaData->east = 180.0;
    metaData->west = -180.0;

    mapProjection = "Undefined";
}

OutFile::~OutFile() {
    free(red);
    free(green);
    free(blue);
    for (size_t i = 0; i < productStuff.size(); i++) {
        delete productStuff[i];
    }
    productStuff.clear();
    free(metaData);
    if (qualityData)
        free(qualityData);
}

void OutFile::setSize(int32_t width, int32_t height) {
    this->width = width;
    this->height = height;
    if (qualityData) {
        free(qualityData);
        qualityData = (uint8_t*) allocateMemory(width, "OutFile::qualityData");
    }
    currentLine = 0;

    for (size_t i = 0; i < productStuff.size(); i++) {
        delete productStuff[i];
    }
    productStuff.clear();
}

int32_t OutFile::getWidth() const {
    return width;
}

int32_t OutFile::getHeight() const {
    return height;
}

void OutFile::setFileName(string fileName) {
    this->fileName = fileName;
}

void OutFile::setPixel(int32_t x, double val, int32_t prod) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile::setPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    productStuff[prod]->lineData[x] = val;
    if (val > fileMaxVal)
        fileMaxVal = val;
    if (val < fileMinVal)
        fileMinVal = val;
}

void OutFile::setPixelRGB(int32_t x, float red, float green, float blue) {
    fprintf(stderr,
            "-E- OutFile::setPixelRGB - RGB not implemented with this file type.\n");
    exit(EXIT_FAILURE);
}

void OutFile::setQuality(int32_t x, uint8_t val) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile::setQuality - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    if (!qualityData) {
        fprintf(stderr, "-E- OutFile::setQuality - qualityData id NULL.\n");
        exit(EXIT_FAILURE);
    }
    qualityData[x] = val;
}

void OutFile::landPixel(int32_t x) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile::landPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    for (size_t prod = 0; prod < productStuff.size(); prod++)
        productStuff[prod]->lineData[x] = landPixelValue;
    if (qualityData)
        qualityData[x] = qualityUnused;
}

void OutFile::fillPixel(int32_t x) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile::fillPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    for (size_t prod = 0; prod < productStuff.size(); prod++)
        productStuff[prod]->lineData[x] = badPixelValue;
    if (qualityData)
        qualityData[x] = qualityUnused;
}

void OutFile::missingPixel(int32_t x) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile::missingPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    for (size_t prod = 0; prod < productStuff.size(); prod++)
        productStuff[prod]->lineData[x] = badPixelValue;
    if (qualityData)
        qualityData[x] = qualityUnused;
}

bool OutFile::setPalette(const char* paletteName, bool applyMask) {
    char* dataRoot;
    string paletteFileName;
    short r[256], g[256], b[256];

    if ((dataRoot = getenv("OCDATAROOT")) == NULL) {
        printf("OCDATAROOT environment variable is not defined.\n");
        return (EXIT_FAILURE);
    }
    paletteFileName = dataRoot;
    paletteFileName += "/common/palette/";
    paletteFileName += paletteName;
    paletteFileName += ".pal";

    if (getlut_file((char*) paletteFileName.c_str(), r, g, b)) {
        fprintf(stderr, "Error reading palette file %s\n",
                paletteFileName.c_str());
        return false;
    }
    if (applyMask) {
        r[landPix] = 160;
        g[landPix] = 82;
        b[landPix] = 45;
        r[fillPix] = 0;
        g[fillPix] = 0;
        b[fillPix] = 0;
    }
    for (int i = 0; i < 256; i++) {
        red[i] = r[i];
        green[i] = g[i];
        blue[i] = b[i];
    }

    return true;
}

void OutFile::setMetaData(meta_l3bType* metaData) {
    *this->metaData = *metaData;
}

/**
 * Add a product for display type output files
 * @param productInfo info structure to copy
 * @return the index for the new product
 */
int32_t OutFile::addProduct(productInfo_t* productInfo) {
    ProductStuff* stuff = new ProductStuff(width, productInfo);

    // setup display scaling
    if (!strcmp(productInfo->displayScale, "linear"))
        stuff->setScale(productInfo->displayMin, productInfo->displayMax,
            Linear);
    else if (!strcmp(productInfo->displayScale, "log"))
        stuff->setScale(productInfo->displayMin, productInfo->displayMax, Log);
    else if (!strcmp(productInfo->displayScale, "arctan"))
        stuff->setScale(productInfo->displayMin, productInfo->displayMax,
            ArcTan);
    else {
        printf("-E- OutFile::addProduct - invalid displayScale = %s\n",
                productInfo->displayScale);
        exit(EXIT_FAILURE);
    }

    productStuff.push_back(stuff);
    return productStuff.size() - 1;
}

int32_t OutFile::addProductNonDisplay(productInfo_t* productInfo) {
    ProductStuff* stuff = new ProductStuff(width, productInfo);

    if (!strcmp(productInfo->dataType, "byte")) {
        stuff->dataStorage = ByteDS;
        stuff->minOutputVal = SCHAR_MIN;
        stuff->maxOutputVal = SCHAR_MAX;
    } else if (!strcmp(productInfo->dataType, "ubyte")) {
        stuff->dataStorage = UByteDS;
        stuff->minOutputVal = 0;
        stuff->maxOutputVal = UCHAR_MAX;
    } else if (!strcmp(productInfo->dataType, "short")) {
        stuff->dataStorage = ShortDS;
        stuff->minOutputVal = SHRT_MIN;
        stuff->maxOutputVal = SHRT_MAX;
    } else if (!strcmp(productInfo->dataType, "ushort")) {
        stuff->dataStorage = UShortDS;
        stuff->minOutputVal = 0;
        stuff->maxOutputVal = USHRT_MAX;
    } else if (!strcmp(productInfo->dataType, "int")) {
        stuff->dataStorage = IntDS;
        stuff->minOutputVal = INT_MIN;
        stuff->maxOutputVal = INT_MAX;
    } else if (!strcmp(productInfo->dataType, "uint")) {
        stuff->dataStorage = UIntDS;
        stuff->minOutputVal = 0;
        stuff->maxOutputVal = UINT_MAX;
    } else if (!strcmp(productInfo->dataType, "float")) {
        stuff->dataStorage = FloatDS;
        stuff->minOutputVal = 0;
        stuff->maxOutputVal = 0;
    } else if (!strcmp(productInfo->dataType, "double")) {
        stuff->dataStorage = DoubleDS;
        stuff->minOutputVal = 0;
        stuff->maxOutputVal = 0;
    } else {
        printf("-E- OutFile::addProductNonDisplay - invalid data type = %s\n",
                productInfo->dataType);
        exit(EXIT_FAILURE);
    }

    // setup scaling
    stuff->setScaleOffset(productInfo->scaleFactor, productInfo->addOffset,
            Linear);
    stuff->missingValue = productInfo->fillValue;

    productStuff.push_back(stuff);
    return productStuff.size() - 1;
}

void OutFile::setMapProjection(string projection) {
    mapProjection = projection;
}

void OutFile::setProj4Info(string projStr, double minX, double maxY) {
    proj4String = projStr;

    for (int i = 0; i < 6; i++)
        tiepoints[i] = 0;
    tiepoints[3] = minX;
    tiepoints[4] = maxY;

    pixscale[0] = resolution;
    pixscale[1] = resolution;
    pixscale[2] = 0;
}

void OutFile::setNumFilledPixels(int32_t num) {
    if (metaData) {
        metaData->data_bins = num;
    }
}

int32_t OutFile::getNumFilledPixels() {
    if (metaData)
        return metaData->data_bins;
    else
        return -1;
}

float OutFile::getPercentFilledPixels() {
    if (metaData) {
        float numPix = width * height;
        return metaData->data_bins / numPix * 100.0;
    } else
        return -1;
}

void OutFile::resetFileMinMax() {
    fileMinVal = DBL_MAX;
    fileMaxVal = 0 - DBL_MAX;
}

void OutFile::setResolution(string resolutionStr) {

    // this will be the # of meters across 1 pixel in the center of the scene
    boost::trim(resolutionStr);
    boost::to_lower(resolutionStr);
    if (resolutionStr.compare("90km") == 0)
        resolution = EARTH_CIRCUMFERENCE / 432.0;
    else if (resolutionStr.compare("36km") == 0)
        resolution = EARTH_CIRCUMFERENCE / 1080.0;
    else if (resolutionStr.compare("18km") == 0)
        resolution = EARTH_CIRCUMFERENCE / 2160.0;
    else if (resolutionStr.compare("9km") == 0)
        resolution = EARTH_CIRCUMFERENCE / 4320.0;
    else if (resolutionStr.compare("4km") == 0)
        resolution = EARTH_CIRCUMFERENCE / 8640.0;
    else if (resolutionStr.compare("2km") == 0)
        resolution = EARTH_CIRCUMFERENCE / 17280.0;
    else if (resolutionStr.compare("1km") == 0)
        resolution = EARTH_CIRCUMFERENCE / 34560.0;
    else if (resolutionStr.compare("hkm") == 0)
        resolution = EARTH_CIRCUMFERENCE / 69120.0;
    else if (resolutionStr.compare("qkm") == 0)
        resolution = EARTH_CIRCUMFERENCE / 138240.0;
    else if (resolutionStr.compare("smi") == 0)
        resolution = EARTH_CIRCUMFERENCE / 4096.0;
    else if (resolutionStr.compare("smi4") == 0)
        resolution = EARTH_CIRCUMFERENCE / 8192.0;
    else if (resolutionStr.compare("land") == 0)
        resolution = EARTH_CIRCUMFERENCE / 8640.0;
    else if (resolutionStr.compare("thirddeg") == 0)
        resolution = EARTH_CIRCUMFERENCE / 1080.0;
    else if (boost::ends_with(resolutionStr, "km")) {
        string val = resolutionStr.substr(0, resolutionStr.length() - 2);
        resolution = atof(val.c_str()) * 1000.0;
    } else if (boost::ends_with(resolutionStr, "deg")) {
        string val = resolutionStr.substr(0, resolutionStr.length() - 3);
        resolution = atof(val.c_str()) / 360.0 * EARTH_CIRCUMFERENCE;
    } else {
        resolution = atof(resolutionStr.c_str());
    }
}

void OutFile::setQualityProcessing(bool val) {
    if (val) {
        // if width is not set yet allocate some dummy memory to flag we want
        // to do SST quality processing.
        if (width <= 0) {
            if (!qualityData)
                qualityData = (uint8_t*) allocateMemory(2,
                    "OutFile::qualityData");
        } else {
            if (qualityData)
                free(qualityData);
            qualityData = (uint8_t*) allocateMemory(width,
                    "OutFile::qualityData");
        }
    } else {
        if (qualityData) {
            free(qualityData);
            qualityData = NULL;
        }
    }
}

bool OutFile::getQualityProcessing() {
    if (qualityData)
        return true;
    else
        return false;
}

//------------------------------------------------------------------------------
// OutFile_pgm
//------------------------------------------------------------------------------

OutFile_pgm::OutFile_pgm() :
OutFile() {
    outfp = NULL;
    fileData = NULL;
}

OutFile_pgm::~OutFile_pgm() {
    if (outfp)
        fclose(outfp);
    if (fileData)
        free(fileData);
}

void OutFile_pgm::setSize(int32_t width, int32_t height) {
    OutFile::setSize(width, height);
    if (fileData)
        free(fileData);
    fileData = (uint8_t*) allocateMemory(width, "OutFile_pgm::fileData");
}

bool OutFile_pgm::open() {
    outfp = fopen(fileName.c_str(), "w");
    if (!outfp)
        return false;

    /* Write pgm header */
    fprintf(outfp, "P5\n");
    fprintf(outfp, "%d %d\n", width, height);
    fprintf(outfp, "255\n");

    return true;
}

void OutFile_pgm::writeLine() {
    for (int i = 0; i < width; i++)
        fileData[i] = (uint8_t) productStuff[0]->calcOutputVal(
            productStuff[0]->lineData[i]);
    fwrite(fileData, 1, width, outfp);
    currentLine++;
}

bool OutFile_pgm::close() {
    fclose(outfp);
    outfp = NULL;
    return true;
}

//------------------------------------------------------------------------------
// OutFile_ppm
//------------------------------------------------------------------------------

void OutFile_ppm::setSize(int32_t width, int32_t height) {
    OutFile::setSize(width, height);
    if (fileData)
        free(fileData);
    fileData = (uint8_t*) allocateMemory(width * 3, "OutFile_ppm::fileData");
}

bool OutFile_ppm::open() {
    outfp = fopen(fileName.c_str(), "w");
    if (!outfp)
        return false;

    /*
     * Write ppm file header
     */
    fprintf(outfp, "P6\n");
    fprintf(outfp, "%d %d\n", width, height);
    fprintf(outfp, "255\n");

    return true;
}

void OutFile_ppm::writeLine() {
    int j = 0;
    for (int i = 0; i < width; i++) {
        uint8_t val = (uint8_t) round(
                productStuff[0]->calcOutputVal(productStuff[0]->lineData[i]));
        fileData[j++] = red[val];
        fileData[j++] = green[val];
        fileData[j++] = blue[val];
    }
    fwrite(fileData, 1, width * 3, outfp);
    currentLine++;
}

//------------------------------------------------------------------------------
// OutFile_ppm_rgb
//------------------------------------------------------------------------------

OutFile_ppm_rgb::OutFile_ppm_rgb() :
OutFile() {
    outfp = NULL;
    fileData = NULL;
}

OutFile_ppm_rgb::~OutFile_ppm_rgb() {
    if (outfp)
        fclose(outfp);
    if (fileData)
        free(fileData);
}

void OutFile_ppm_rgb::setSize(int32_t width, int32_t height) {
    OutFile::setSize(width, height);
    if (fileData)
        free(fileData);
    fileData = (uint8_t*) allocateMemory(width * 3,
            "OutFile_ppm_rgb::fileData");
}

bool OutFile_ppm_rgb::open() {
    outfp = fopen(fileName.c_str(), "w");
    if (!outfp)
        return false;

    /*
     * Write ppm file header
     */
    fprintf(outfp, "P6\n");
    fprintf(outfp, "%d %d\n", width, height);
    fprintf(outfp, "255\n");
    return true;
}

bool OutFile_ppm_rgb::close() {
    fclose(outfp);
    outfp = NULL;
    return true;
}

void OutFile_ppm_rgb::setPixel(int32_t x, double val, int32_t prod) {
    fprintf(stderr,
            "-E- OutFile_ppm_rgb::setPixel - only RGB is implemented with this file type.\n");
    exit(EXIT_FAILURE);
}

void OutFile_ppm_rgb::setPixelRGB(int32_t x, float red, float green,
        float blue) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile_ppm_rgb::setPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }

    uint8_t *ptr = fileData + x * 3;

    *ptr = round(productStuff[0]->calcOutputVal(red));
    ptr++;
    *ptr = round(productStuff[0]->calcOutputVal(green));
    ptr++;
    *ptr = round(productStuff[0]->calcOutputVal(blue));

    // do this to keep the file min/max reasonable
    if (red > fileMaxVal)
        fileMaxVal = red;
    if (red < fileMinVal)
        fileMinVal = red;
}

void OutFile_ppm_rgb::landPixel(int32_t x) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile_ppm_rgb::landPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    uint8_t *ptr = fileData + x * 3;
    *ptr = landPix;
    ptr++;
    *ptr = landPix;
    ptr++;
    *ptr = landPix;
}

void OutFile_ppm_rgb::fillPixel(int32_t x) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile_ppm_rgb::fillPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    uint8_t *ptr = fileData + x * 3;
    *ptr = fillPix;
    ptr++;
    *ptr = fillPix;
    ptr++;
    *ptr = fillPix;
}

void OutFile_ppm_rgb::missingPixel(int32_t x) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile_ppm_rgb::missingPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    uint8_t *ptr = fileData + x * 3;
    *ptr = fillPix;
    ptr++;
    *ptr = fillPix;
    ptr++;
    *ptr = fillPix;
}

void OutFile_ppm_rgb::writeLine() {
    fwrite(fileData, 1, width * 3, outfp);
    currentLine++;
}

//------------------------------------------------------------------------------
// OutFile_png
//------------------------------------------------------------------------------

OutFile_png::OutFile_png(bool color) :
OutFile() {
    isColor = color;
    outfp = NULL;
    fileData = NULL;
    info_ptr = NULL;
    png_ptr = NULL;
}

OutFile_png::~OutFile_png() {
    if (outfp)
        fclose(outfp);
    if (fileData)
        free(fileData);
}

void OutFile_png::setSize(int32_t width, int32_t height) {
    OutFile::setSize(width, height);
    if (fileData)
        free(fileData);
    fileData = (uint8_t*) allocateMemory(width, "OutFile_png::fileData");
}

bool OutFile_png::open() {
    outfp = fopen(fileName.c_str(), "w");
    if (!outfp)
        return false;

    png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING,
            NULL, NULL, NULL);
    if (!png_ptr) {
        fprintf(stderr, "-E- Unable to create PNG write structure.\n");
        exit(EXIT_FAILURE);
    }

    info_ptr = png_create_info_struct(png_ptr);
    if (!info_ptr) {
        png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
        fprintf(stderr, "-E- Unable to create PNG info structure.\n");
        exit(EXIT_FAILURE);
    }
    if (setjmp(png_jmpbuf(png_ptr))) {
        png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
        fprintf(stderr, "-E- Unable to call PNG setjmp().\n");
        exit(EXIT_FAILURE);
    }
    png_init_io(png_ptr, outfp);

    if (isColor) {

        // color palette
        png_set_IHDR(png_ptr, info_ptr, width, height, 8,
                PNG_COLOR_TYPE_PALETTE, PNG_INTERLACE_NONE,
                PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);

        uint8_t pal[256 * 3];
        for (int i = 0; i < 256; i++) {
            pal[i * 3] = red[i];
            pal[i * 3 + 1] = green[i];
            pal[i * 3 + 2] = blue[i];
        }
        png_set_PLTE(png_ptr, info_ptr, (png_const_colorp) pal, 256);
    } else {

        // Grayscale
        png_set_IHDR(png_ptr, info_ptr, width, height, 8, PNG_COLOR_TYPE_GRAY,
                PNG_INTERLACE_NONE,
                PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
    }
    png_write_info(png_ptr, info_ptr);

    return true;
}

void OutFile_png::writeLine() {
    for (int i = 0; i < width; i++)
        fileData[i] = (uint8_t) round(
            productStuff[0]->calcOutputVal(productStuff[0]->lineData[i]));
    png_write_row(png_ptr, (png_bytep) fileData);
    currentLine++;
}

bool OutFile_png::close() {
    png_write_end(png_ptr, info_ptr);
    png_destroy_write_struct(&png_ptr, &info_ptr);

    fclose(outfp);
    outfp = NULL;
    return true;
}

//------------------------------------------------------------------------------
// OutFile_png_rgb
//------------------------------------------------------------------------------

OutFile_png_rgb::OutFile_png_rgb() :
OutFile() {
    outfp = NULL;
    fileData = NULL;
    info_ptr = NULL;
    png_ptr = NULL;
}

OutFile_png_rgb::~OutFile_png_rgb() {
    if (outfp)
        fclose(outfp);
    if (fileData)
        free(fileData);
}

void OutFile_png_rgb::setSize(int32_t width, int32_t height) {
    OutFile::setSize(width, height);
    if (fileData)
        free(fileData);
    fileData = (uint8_t*) allocateMemory(width * 3,
            "OutFile_png_rgb::fileData");
}

bool OutFile_png_rgb::open() {
    outfp = fopen(fileName.c_str(), "w");
    if (!outfp)
        return false;

    png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING,
            NULL, NULL, NULL);
    if (!png_ptr) {
        fprintf(stderr, "-E- Unable to create PNG write structure.\n");
        exit(EXIT_FAILURE);
    }

    info_ptr = png_create_info_struct(png_ptr);
    if (!info_ptr) {
        png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
        fprintf(stderr, "-E- Unable to create PNG info structure.\n");
        exit(EXIT_FAILURE);
    }
    if (setjmp(png_jmpbuf(png_ptr))) {
        png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
        fprintf(stderr, "-E- Unable to call PNG setjmp().\n");
        exit(EXIT_FAILURE);
    }
    png_init_io(png_ptr, outfp);
    png_set_IHDR(png_ptr, info_ptr, width, height, 8, PNG_COLOR_TYPE_RGB,
            PNG_INTERLACE_NONE,
            PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
    png_write_info(png_ptr, info_ptr);

    return true;
}

bool OutFile_png_rgb::close() {
    if (setjmp(png_jmpbuf(png_ptr))) {
        fprintf(stderr,
                "-E- OutFile_png_rgb::close - Unable to call PNG setjmp().\n");
        exit(EXIT_FAILURE);
    }
    png_write_end(png_ptr, info_ptr);
    png_destroy_write_struct(&png_ptr, &info_ptr);

    fclose(outfp);
    outfp = NULL;
    return true;
}

void OutFile_png_rgb::setPixel(int32_t x, double val, int32_t prod) {
    fprintf(stderr,
            "-E- OutFile_png_rgb::setPixel - only RGB is implemented with this file type.\n");
    exit(EXIT_FAILURE);
}

void OutFile_png_rgb::setPixelRGB(int32_t x, float red, float green,
        float blue) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile_png_rgb::setPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }

    uint8_t *ptr = fileData + x * 3;

    *ptr = round(productStuff[0]->calcOutputVal(red));
    ptr++;
    *ptr = round(productStuff[0]->calcOutputVal(green));
    ptr++;
    *ptr = round(productStuff[0]->calcOutputVal(blue));

    // do this to keep the file min/max reasonable
    if (red > fileMaxVal)
        fileMaxVal = red;
    if (red < fileMinVal)
        fileMinVal = red;
}

void OutFile_png_rgb::landPixel(int32_t x) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile_png_rgb::landPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    uint8_t *ptr = fileData + x * 3;
    *ptr = landPix;
    ptr++;
    *ptr = landPix;
    ptr++;
    *ptr = landPix;
}

void OutFile_png_rgb::fillPixel(int32_t x) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile_png_rgb::fillPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    uint8_t *ptr = fileData + x * 3;
    *ptr = fillPix;
    ptr++;
    *ptr = fillPix;
    ptr++;
    *ptr = fillPix;
}

void OutFile_png_rgb::missingPixel(int32_t x) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile_png_rgb::missingPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    uint8_t *ptr = fileData + x * 3;
    *ptr = fillPix;
    ptr++;
    *ptr = fillPix;
    ptr++;
    *ptr = fillPix;
}

void OutFile_png_rgb::writeLine() {

    if (setjmp(png_jmpbuf(png_ptr))) {
        fprintf(stderr,
                "-E- OutFile_png_rgb::writeLine - Unable to call PNG setjmp().\n");
        exit(EXIT_FAILURE);
    }

    png_write_row(png_ptr, (png_bytep) fileData);
    currentLine++;
}

//------------------------------------------------------------------------------
// OutFile_tiff
//------------------------------------------------------------------------------

OutFile_tiff::~OutFile_tiff() {
    close();
}

bool OutFile_tiff::open() {
    currentLine = 0;

    // open TIFF file
    tiff = XTIFFOpen(fileName.c_str(), "w");
    if (tiff == NULL) {
        cerr << "-E- Could not open outgoing TIFF image" << endl;
        exit(EXIT_FAILURE);
    }
    TIFFSetField(tiff, TIFFTAG_IMAGEWIDTH, width);
    TIFFSetField(tiff, TIFFTAG_IMAGELENGTH, height);
    TIFFSetField(tiff, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
    TIFFSetField(tiff, TIFFTAG_COMPRESSION, COMPRESSION_LZW);
    TIFFSetField(tiff, TIFFTAG_PREDICTOR, PREDICTOR_NONE);
    TIFFSetField(tiff, TIFFTAG_ROWSPERSTRIP, height);

    // open GeoTIFF interface
    bool hasGeotiffInfo = false;
    gtif = GTIFNew(tiff);
    if (gtif == NULL) {
        cerr << "-E- Could not create geoTIFF structure" << endl;
        exit(EXIT_FAILURE);
    }

    // define GeoTIFF keys for lat/lon projection
    if (mapProjection.compare("Equidistant Cylindrical") == 0) {

        double north, south, east, west;
        if (metaData) {
            north = metaData->north;
            south = metaData->south;
            east = metaData->east;
            west = metaData->west;
        } else {
            north = 90.0;
            south = -90.0;
            east = 180.0;
            west = -180.0;
        }

        // pixel width, height in degrees
        pixscale[0] = (east - west) / width;
        pixscale[1] = (north - south) / height;

        // top left corner pixel lat, lon
        for (int i = 0; i < 6; i++)
            tiepoints[i] = 0;
        tiepoints[3] = west;
        tiepoints[4] = north;

        GTIFKeySet(gtif, GTModelTypeGeoKey, TYPE_SHORT, 1, ModelGeographic);
        GTIFKeySet(gtif, GeographicTypeGeoKey, TYPE_SHORT, 1, GCS_WGS_84);
        hasGeotiffInfo = true;
    }
    // otherwise, parse the proj4 string
    else
        hasGeotiffInfo = (GTIFSetFromProj4(gtif, proj4String.c_str()));

    if (hasGeotiffInfo) {

        // write GeoTIFF keys
        GTIFKeySet(gtif, GTRasterTypeGeoKey, TYPE_SHORT, 1, RasterPixelIsArea);
        GTIFWriteKeys(gtif);

        // write GeoTIFF tags in map units
        TIFFSetField(tiff, GTIFF_PIXELSCALE, 3, pixscale);
        TIFFSetField(tiff, GTIFF_TIEPOINTS, 6, tiepoints);
    }

    // define colormap
    setTiffColor();
    return true;
}

bool OutFile_tiff::close() {
    if (gtif) {
        GTIFFree(gtif);
        XTIFFClose(tiff);
        tiff = NULL;
        gtif = NULL;
    }
    return true;
}

//----- OutFile_tiff_color -----

OutFile_tiff_color::~OutFile_tiff_color() {
    if (fileData)
        free(fileData);
}

void OutFile_tiff_color::setSize(int32_t width, int32_t height) {
    OutFile::setSize(width, height);
    if (fileData)
        free(fileData);
    fileData = (uint8_t*) allocateMemory(width, "OutFile_tiff_color::fileData");
}

void OutFile_tiff_color::writeLine() {
    for (int i = 0; i < width; i++)
        fileData[i] = (uint8_t) round(
                productStuff[0]->calcOutputVal(productStuff[0]->lineData[i]));

    if (TIFFWriteScanline(tiff, (void*) fileData, currentLine) == 0) {
        cerr << "-E- Could not write TIFF image line " << currentLine << endl;
        exit(EXIT_FAILURE);
    }
    currentLine++;
}

void OutFile_tiff_color::setTiffColor() {

    TIFFSetField(tiff, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_PALETTE);
    TIFFSetField(tiff, TIFFTAG_BITSPERSAMPLE, 8);
    TIFFSetField(tiff, TIFFTAG_SAMPLESPERPIXEL, 1);

    int ncolors = 256; // convert bytes to short
    uint16_t* rr = (uint16_t*) malloc(ncolors * sizeof(uint16_t));
    uint16_t* gg = (uint16_t*) malloc(ncolors * sizeof(uint16_t));
    uint16_t* bb = (uint16_t*) malloc(ncolors * sizeof(uint16_t));
    if (rr == NULL || gg == NULL || bb == NULL) {
        cerr << "-E- Could not allocate memory for TIFF color map" << endl;
        exit(EXIT_FAILURE);
    }
    for (int i = 0; i < ncolors; i++) {
        rr[i] = red[i] << 8;
        gg[i] = green[i] << 8;
        bb[i] = blue[i] << 8;
    }
    TIFFSetField(tiff, TIFFTAG_COLORMAP, rr, gg, bb);
    free(rr);
    free(gg);
    free(bb);

}

//----- OutFile_tiff_gray -----

OutFile_tiff_gray::~OutFile_tiff_gray() {
    if (fileData)
        free(fileData);
}

void OutFile_tiff_gray::setSize(int32_t width, int32_t height) {
    OutFile::setSize(width, height);
    if (fileData)
        free(fileData);
    fileData = (float*) allocateMemory(width * sizeof(float),
            "OutFile_tiff_gray::fileData");
}

void OutFile_tiff_gray::writeLine() {
    for (int i = 0; i < width; i++)
        fileData[i] = productStuff[0]->lineData[i];

    if (TIFFWriteScanline(tiff, (void*) fileData, currentLine) == 0) {
        cerr << "-E- Could not write TIFF image line " << currentLine << endl;
        exit(EXIT_FAILURE);
    }
    currentLine++;
}

void OutFile_tiff_gray::setTiffColor() {
    TIFFSetField(tiff, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK);
    TIFFSetField(tiff, TIFFTAG_BITSPERSAMPLE, 32);
    TIFFSetField(tiff, TIFFTAG_SAMPLESPERPIXEL, 1);
    TIFFSetField(tiff, TIFFTAG_SAMPLEFORMAT, SAMPLEFORMAT_IEEEFP);
}

//----- OutFile_tiff_rgb -----

OutFile_tiff_rgb::~OutFile_tiff_rgb() {
    if (fileData)
        free(fileData);
}

void OutFile_tiff_rgb::setSize(int32_t width, int32_t height) {
    OutFile::setSize(width, height);
    if (fileData)
        free(fileData);
    fileData = (uint8_t*) allocateMemory(width * 3,
            "OutFile_tiff_rgb::fileData");
}

void OutFile_tiff_rgb::writeLine() {
    if (TIFFWriteScanline(tiff, (void*) fileData, currentLine) == 0) {
        cerr << "-E- Could not write TIFF image line " << currentLine << endl;
        exit(EXIT_FAILURE);
    }
    currentLine++;
}

void OutFile_tiff_rgb::setTiffColor() {
    TIFFSetField(tiff, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
    TIFFSetField(tiff, TIFFTAG_BITSPERSAMPLE, 8);
    TIFFSetField(tiff, TIFFTAG_SAMPLESPERPIXEL, 3);
    TIFFSetField(tiff, TIFFTAG_COMPRESSION, COMPRESSION_DEFLATE);
}

void OutFile_tiff_rgb::setPixel(int32_t x, double val, int32_t prod) {
    fprintf(stderr,
            "-E- OutFile_tiff_rgb::setPixel - only RGB is implemented with this file type.\n");
    exit(EXIT_FAILURE);
}

void OutFile_tiff_rgb::setPixelRGB(int32_t x, float red, float green,
                                   float blue) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile_tiff_rgb::setPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }

    uint8_t *ptr = fileData + x * 3;

    *ptr = round(productStuff[0]->calcOutputVal(red));
    ptr++;
    *ptr = round(productStuff[0]->calcOutputVal(green));
    ptr++;
    *ptr = round(productStuff[0]->calcOutputVal(blue));

    if (red > fileMaxVal) // keep file min/max reasonable
        fileMaxVal = red;
    if (red < fileMinVal)
        fileMinVal = red;
}

void OutFile_tiff_rgb::landPixel(int32_t x) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile_tiff_rgb::landPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    uint8_t *ptr = fileData + x * 3;
    *ptr = landPix;
    ptr++;
    *ptr = landPix;
    ptr++;
    *ptr = landPix;
}

void OutFile_tiff_rgb::fillPixel(int32_t x) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile_tiff_rgb::fillPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    uint8_t *ptr = fileData + x * 3;
    *ptr = fillPix;
    ptr++;
    *ptr = fillPix;
    ptr++;
    *ptr = fillPix;
}

void OutFile_tiff_rgb::missingPixel(int32_t x) {
    if (x < 0 || x >= width) {
        fprintf(stderr,
                "-E- OutFile_tiff_rgb::missingPixel - x=%d is not within range, width=%d.\n",
                x, width);
        exit(EXIT_FAILURE);
    }
    uint8_t *ptr = fileData + x * 3;
    *ptr = fillPix;
    ptr++;
    *ptr = fillPix;
    ptr++;
    *ptr = fillPix;
}

//------------------------------------------------------------------------------
// OutFile_hdf4
//------------------------------------------------------------------------------

OutFile_hdf4::OutFile_hdf4() :
OutFile() {
    fileData = NULL;
    sdfid = -1;
    sdsid = -1;
    quality_sdsid = -1;
    hdfDataType = DFNT_FLOAT32;
}

OutFile_hdf4::~OutFile_hdf4() {
    if (fileData)
        free(fileData);
}

void OutFile_hdf4::setSize(int32_t width, int32_t height) {
    OutFile::setSize(width, height);
    if (fileData)
        free(fileData);
    fileData = NULL;
}

bool OutFile_hdf4::open() {
    const char* tmpStr;
    float tmpFloat;
    int32_t tmpInt32;

    currentLine = 0;

    sdfid = SDstart(fileName.c_str(), DFACC_CREATE);

    if (sdfid < 0) {
        printf("-E- Could not create HDF4 file %s\n", fileName.c_str());
        exit(EXIT_FAILURE);
    }

    get_time(metaData->ptime);

    string prodName;
    size_t pos = fileName.find_last_of('/');
    if (pos == string::npos)
        prodName = fileName;
    else
        prodName = fileName.substr(pos + 1);
    DPTB(
            SDsetattr(sdfid, "Product Name", DFNT_CHAR, prodName.size() + 1, (VOIDP) prodName.c_str()));
    DPTB(
            SDsetattr(sdfid, "Sensor Name", DFNT_CHAR, strlen(metaData->sensor_name) + 1, (VOIDP) metaData->sensor_name));
    DPTB(
            SDsetattr(sdfid, "Sensor", DFNT_CHAR, strlen(metaData->sensor) + 1, (VOIDP) metaData->sensor));
    DPTB(
            SDsetattr(sdfid, "Title", DFNT_CHAR, strlen(metaData->title) + 1, (VOIDP) metaData->title));
    DPTB(
            SDsetattr(sdfid, "Data Center", DFNT_CHAR, strlen(metaData->data_center) + 1, (VOIDP) metaData->data_center));
    DPTB(
            SDsetattr(sdfid, "Station Name", DFNT_CHAR, strlen(metaData->station) + 1, (VOIDP) metaData->station));
    DPTB(
            SDsetattr(sdfid, "Station Latitude", DFNT_FLOAT32, 1, (VOIDP) & metaData->station_lat));
    DPTB(
            SDsetattr(sdfid, "Station Longitude", DFNT_FLOAT32, 1, (VOIDP) & metaData->station_lon));
    DPTB(
            SDsetattr(sdfid, "Mission", DFNT_CHAR, strlen(metaData->mission) + 1, (VOIDP) metaData->mission));
    DPTB(
            SDsetattr(sdfid, "Mission Characteristics", DFNT_CHAR, strlen(metaData->mission_char) + 1, (VOIDP) metaData->mission_char));
    DPTB(
            SDsetattr(sdfid, "Sensor Characteristics", DFNT_CHAR, strlen(metaData->sensor_char) + 1, (VOIDP) metaData->sensor_char));
    DPTB(
            SDsetattr(sdfid, "Product Type", DFNT_CHAR, strlen(metaData->prod_type) + 1, (VOIDP) metaData->prod_type));
    DPTB(
            SDsetattr(sdfid, "Processing Version", DFNT_CHAR, strlen(metaData->pversion) + 1, (VOIDP) metaData->pversion));
    DPTB(
            SDsetattr(sdfid, "Software Name", DFNT_CHAR, strlen(metaData->soft_name) + 1, (VOIDP) metaData->soft_name));
    DPTB(
            SDsetattr(sdfid, "Software Version", DFNT_CHAR, strlen(metaData->soft_ver) + 1, (VOIDP) metaData->soft_ver));
    DPTB(
            SDsetattr(sdfid, "Processing Time", DFNT_CHAR, strlen(metaData->ptime) + 1, (VOIDP) metaData->ptime));
    DPTB(
            SDsetattr(sdfid, "Input Files", DFNT_CHAR, strlen(metaData->infiles) + 1, (VOIDP) metaData->infiles));
    DPTB(
            SDsetattr(sdfid, "Processing Control", DFNT_CHAR, strlen(metaData->proc_con) + 1, (VOIDP) metaData->proc_con));
    DPTB(
            SDsetattr(sdfid, "Input Parameters", DFNT_CHAR, strlen(metaData->input_parms) + 1, (VOIDP) metaData->input_parms));
    DPTB(
            SDsetattr(sdfid, "L2 Flag Names", DFNT_CHAR, strlen(metaData->flag_names) + 1, (VOIDP) metaData->flag_names));

    short syear, sday, eyear, eday;
    double ssec, esec;
    int32_t smsec, emsec;
    unix2yds(metaData->startTime, &syear, &sday, &ssec);
    smsec = (int32_t) (ssec * 1000.0);
    unix2yds(metaData->endTime, &eyear, &eday, &esec);
    emsec = (int32_t) (esec * 1000.0);
    DPTB(SDsetattr(sdfid, "Period Start Year", DFNT_INT16, 1, (VOIDP) & syear));
    DPTB(SDsetattr(sdfid, "Period Start Day", DFNT_INT16, 1, (VOIDP) & sday));
    DPTB(SDsetattr(sdfid, "Period End Year", DFNT_INT16, 1, (VOIDP) & eyear));
    DPTB(SDsetattr(sdfid, "Period End Day", DFNT_INT16, 1, (VOIDP) & eday));
    tmpStr = ydhmsf(metaData->startTime, 'G');
    DPTB(
            SDsetattr(sdfid, "Start Time", DFNT_CHAR, strlen(tmpStr) + 1, (VOIDP) tmpStr));
    tmpStr = ydhmsf(metaData->endTime, 'G');
    DPTB(
            SDsetattr(sdfid, "End Time", DFNT_CHAR, strlen(tmpStr) + 1, (VOIDP) tmpStr));
    DPTB(SDsetattr(sdfid, "Start Year", DFNT_INT16, 1, (VOIDP) & syear));
    DPTB(SDsetattr(sdfid, "Start Day", DFNT_INT16, 1, (VOIDP) & sday));
    DPTB(SDsetattr(sdfid, "Start Millisec", DFNT_INT32, 1, (VOIDP) & smsec));
    DPTB(SDsetattr(sdfid, "End Year", DFNT_INT16, 1, (VOIDP) & eyear));
    DPTB(SDsetattr(sdfid, "End Day", DFNT_INT16, 1, (VOIDP) & eday));
    DPTB(SDsetattr(sdfid, "End Millisec", DFNT_INT32, 1, (VOIDP) & emsec));

    DPTB(
            SDsetattr(sdfid, "Start Orbit", DFNT_INT32, 1, (VOIDP) & metaData->start_orb));
    DPTB(
            SDsetattr(sdfid, "End Orbit", DFNT_INT32, 1, (VOIDP) & metaData->end_orb));
    DPTB(SDsetattr(sdfid, "Orbit", DFNT_INT32, 1, (VOIDP) & metaData->orbit));
    DPTB(
            SDsetattr(sdfid, "Map Projection", DFNT_CHAR, mapProjection.size() + 1, mapProjection.c_str()));
    DPTB(
            SDsetattr(sdfid, "Latitude Units", DFNT_CHAR, 14, (VOIDP) "degrees North"));
    DPTB(
            SDsetattr(sdfid, "Longitude Units", DFNT_CHAR, 13, (VOIDP) "degrees East"));
    DPTB(
            SDsetattr(sdfid, "Northernmost Latitude", DFNT_FLOAT32, 1, (VOIDP) & metaData->north));
    DPTB(
            SDsetattr(sdfid, "Southernmost Latitude", DFNT_FLOAT32, 1, (VOIDP) & metaData->south));
    DPTB(
            SDsetattr(sdfid, "Westernmost Longitude", DFNT_FLOAT32, 1, (VOIDP) & metaData->west));
    DPTB(
            SDsetattr(sdfid, "Easternmost Longitude", DFNT_FLOAT32, 1, (VOIDP) & metaData->east));

    float latStep = (metaData->north - metaData->south) / (float) getHeight();
    DPTB(SDsetattr(sdfid, "Latitude Step", DFNT_FLOAT32, 1, (VOIDP) & latStep));
    float lonStep;
    if (metaData->east < metaData->west)
        lonStep = (360 + metaData->east - metaData->west) / (float) getWidth();
    else
        lonStep = (metaData->east - metaData->west) / (float) getWidth();
    DPTB(SDsetattr(sdfid, "Longitude Step", DFNT_FLOAT32, 1, (VOIDP) & lonStep));

    tmpFloat = metaData->south + latStep / 2.0;
    DPTB(
            SDsetattr(sdfid, "SW Point Latitude", DFNT_FLOAT32, 1, (VOIDP) & tmpFloat));
    tmpFloat = metaData->west + lonStep / 2.0;
    DPTB(
            SDsetattr(sdfid, "SW Point Longitude", DFNT_FLOAT32, 1, (VOIDP) & tmpFloat));

    DPTB(
            SDsetattr(sdfid, "Data Bins", DFNT_INT32, 1, (VOIDP) & metaData->data_bins));

    tmpInt32 = getHeight();
    DPTB(SDsetattr(sdfid, "Number of Lines", DFNT_INT32, 1, (VOIDP) & tmpInt32));
    tmpInt32 = getWidth();
    DPTB(
            SDsetattr(sdfid, "Number of Columns", DFNT_INT32, 1, (VOIDP) & tmpInt32));
    DPTB(
            SDsetattr(sdfid, "Parameter", DFNT_CHAR, strlen(productStuff[0]->productInfo->description) + 1, (VOIDP) productStuff[0]->productInfo->description));
    DPTB(SDsetattr(sdfid, "Measure", DFNT_CHAR, 5, (VOIDP) "Mean"));
    DPTB(
            SDsetattr(sdfid, "Units", DFNT_CHAR, strlen(productStuff[0]->productInfo->units) + 1, (VOIDP) productStuff[0]->productInfo->units));

    // we only use linear scaling for data storage
    tmpStr = "linear";
    DPTB(
            SDsetattr(sdfid, "Scaling", DFNT_CHAR, strlen(tmpStr) + 1, (VOIDP) tmpStr));
    tmpStr = "(Slope*l3m_data) + Intercept = Parameter value";
    DPTB(
            SDsetattr(sdfid, "Scaling Equation", DFNT_CHAR, strlen(tmpStr) + 1, (VOIDP) tmpStr));

    float tmpScale;
    float tmpOffset;
    const char* imageScalingApplied;

    switch (productStuff[0]->dataStorage) {
    case FloatDS:
    case DoubleDS:
        tmpScale = 1.0;
        tmpOffset = 0.0;
        imageScalingApplied = "No";
        break;
    default:
        tmpScale = productStuff[0]->scale;
        tmpOffset = productStuff[0]->offset;
        if (tmpScale == 1.0 && tmpOffset == 0.0)
            imageScalingApplied = "No";
        else
            imageScalingApplied = "Yes";
    }

    DPTB(SDsetattr(sdfid, "Slope", DFNT_FLOAT32, 1, (VOIDP) & tmpScale));
    DPTB(SDsetattr(sdfid, "Intercept", DFNT_FLOAT32, 1, (VOIDP) & tmpOffset));

    tmpFloat = productStuff[0]->productInfo->validMin;
    DPTB(SDsetattr(sdfid, "Data Minimum", DFNT_FLOAT32, 1, (VOIDP) & tmpFloat));
    tmpFloat = productStuff[0]->productInfo->validMax;
    DPTB(SDsetattr(sdfid, "Data Maximum", DFNT_FLOAT32, 1, (VOIDP) & tmpFloat));
    tmpFloat = productStuff[0]->productInfo->displayMin;
    DPTB(
            SDsetattr(sdfid, "Suggested Image Scaling Minimum", DFNT_FLOAT32, 1, (VOIDP) & tmpFloat));
    tmpFloat = productStuff[0]->productInfo->displayMax;
    DPTB(
            SDsetattr(sdfid, "Suggested Image Scaling Maximum", DFNT_FLOAT32, 1, (VOIDP) & tmpFloat));

    tmpStr = strdup(productStuff[0]->productInfo->displayScale);
    upcase((char*) tmpStr);
    DPTB(
            SDsetattr(sdfid, "Suggested Image Scaling Type", DFNT_CHAR, strlen(tmpStr) + 1, (VOIDP) tmpStr));
    free((void*) tmpStr);

    DPTB(
            SDsetattr(sdfid, "Suggested Image Scaling Applied", DFNT_CHAR, strlen(imageScalingApplied) + 1, (VOIDP) imageScalingApplied));
    DPTB(
            SDsetattr(sdfid, "_lastModified", DFNT_CHAR, strlen(metaData->ptime) + 1, (VOIDP) metaData->ptime));

    // delete file data
    if (fileData)
        free(fileData);

    int32_t dims[2];
    dims[0] = height;
    dims[1] = width;

    if (!strcmp(productStuff[0]->productInfo->dataType, "byte")) {
        hdfDataType = DFNT_INT8;
        sdsid = SDcreate(sdfid, "l3m_data", hdfDataType, 2, dims);
        int8_t tmp = productStuff[0]->missingValue;
        DPTB(SDsetattr(sdsid, "Fill", hdfDataType, 1, (VOIDP) & tmp));
        fileData = allocateMemory(width * sizeof (int8_t),
                "OutFile_hdf4::open fileData");
    } else if (!strcmp(productStuff[0]->productInfo->dataType, "ubyte")) {
        hdfDataType = DFNT_UINT8;
        sdsid = SDcreate(sdfid, "l3m_data", hdfDataType, 2, dims);
        uint8_t tmp = productStuff[0]->missingValue;
        DPTB(SDsetattr(sdsid, "Fill", hdfDataType, 1, (VOIDP) & tmp));
        fileData = allocateMemory(width * sizeof (uint8_t),
                "OutFile_hdf4::open fileData");
    } else if (!strcmp(productStuff[0]->productInfo->dataType, "short")) {
        hdfDataType = DFNT_INT16;
        sdsid = SDcreate(sdfid, "l3m_data", hdfDataType, 2, dims);
        int16_t tmp = productStuff[0]->missingValue;
        DPTB(SDsetattr(sdsid, "Fill", hdfDataType, 1, (VOIDP) & tmp));
        fileData = allocateMemory(width * sizeof (int16_t),
                "OutFile_hdf4::open fileData");
    } else if (!strcmp(productStuff[0]->productInfo->dataType, "ushort")) {
        hdfDataType = DFNT_UINT16;
        sdsid = SDcreate(sdfid, "l3m_data", hdfDataType, 2, dims);
        uint16_t tmp = productStuff[0]->missingValue;
        DPTB(SDsetattr(sdsid, "Fill", hdfDataType, 1, (VOIDP) & tmp));
        fileData = allocateMemory(width * sizeof (uint16_t),
                "OutFile_hdf4::open fileData");
    } else if (!strcmp(productStuff[0]->productInfo->dataType, "int")) {
        hdfDataType = DFNT_INT32;
        sdsid = SDcreate(sdfid, "l3m_data", hdfDataType, 2, dims);
        int32_t tmp = productStuff[0]->missingValue;
        DPTB(SDsetattr(sdsid, "Fill", hdfDataType, 1, (VOIDP) & tmp));
        fileData = allocateMemory(width * sizeof (int32_t),
                "OutFile_hdf4::open fileData");
    } else if (!strcmp(productStuff[0]->productInfo->dataType, "uint")) {
        hdfDataType = DFNT_UINT32;
        sdsid = SDcreate(sdfid, "l3m_data", hdfDataType, 2, dims);
        uint32_t tmp = productStuff[0]->missingValue;
        DPTB(SDsetattr(sdsid, "Fill", hdfDataType, 1, (VOIDP) & tmp));
        fileData = allocateMemory(width * sizeof (uint32_t),
                "OutFile_hdf4::open fileData");
    } else if (!strcmp(productStuff[0]->productInfo->dataType, "float")) {
        hdfDataType = DFNT_FLOAT32;
        sdsid = SDcreate(sdfid, "l3m_data", hdfDataType, 2, dims);
        float tmp = productStuff[0]->missingValue;
        DPTB(SDsetattr(sdsid, "Fill", hdfDataType, 1, (VOIDP) & tmp));
        fileData = allocateMemory(width * sizeof (float),
                "OutFile_hdf4::open fileData");
    } else if (!strcmp(productStuff[0]->productInfo->dataType, "double")) {
        hdfDataType = DFNT_FLOAT64;
        sdsid = SDcreate(sdfid, "l3m_data", hdfDataType, 2, dims);
        double tmp = productStuff[0]->missingValue;
        DPTB(SDsetattr(sdsid, "Fill", hdfDataType, 1, (VOIDP) & tmp));
        fileData = allocateMemory(width * sizeof (double),
                "OutFile_hdf4::open fileData");
    } else {
        printf("-E- Data type %s, not supported\n",
                productStuff[0]->productInfo->dataType);
        exit(EXIT_FAILURE);
    }

    tmpStr = "linear";
    DPTB(
            SDsetattr(sdsid, "Scaling", DFNT_CHAR, strlen(tmpStr) + 1, (VOIDP) tmpStr));
    tmpStr = "(Slope*l3m_data) + Intercept = Parameter value";
    DPTB(
            SDsetattr(sdsid, "Scaling Equation", DFNT_CHAR, strlen(tmpStr) + 1, (VOIDP) tmpStr));
    DPTB(SDsetattr(sdsid, "Slope", DFNT_FLOAT32, 1, (VOIDP) & tmpScale));
    DPTB(SDsetattr(sdsid, "Intercept", DFNT_FLOAT32, 1, (VOIDP) & tmpOffset));

    // create the SST quality data set
    if (qualityData) {
        quality_sdsid = SDcreate(sdfid, "l3m_qual", DFNT_UINT8, 2, dims);
        int32_t validRange[2];
        validRange[0] = 0;
        validRange[1] = 2;
        DPTB(
                SDsetattr(quality_sdsid, "valid_range", DFNT_INT32, 2, (VOIDP) validRange));
    }

    return true;
}

void OutFile_hdf4::writeLine() {

    switch (productStuff[0]->dataStorage) {
    case ByteDS:
        for (int i = 0; i < width; i++)
            ((int8_t*) fileData)[i] = round(
                productStuff[0]->calcOutputVal(
                productStuff[0]->lineData[i]));
        break;
    case UByteDS:
        for (int i = 0; i < width; i++)
            ((uint8_t*) fileData)[i] = round(
                productStuff[0]->calcOutputVal(
                productStuff[0]->lineData[i]));
        break;
    case ShortDS:
        for (int i = 0; i < width; i++)
            ((int16_t*) fileData)[i] = round(
                productStuff[0]->calcOutputVal(
                productStuff[0]->lineData[i]));
        break;
    case UShortDS:
        for (int i = 0; i < width; i++)
            ((uint16_t*) fileData)[i] = round(
                productStuff[0]->calcOutputVal(
                productStuff[0]->lineData[i]));
        break;
    case IntDS:
        for (int i = 0; i < width; i++)
            ((int32_t*) fileData)[i] = round(
                productStuff[0]->calcOutputVal(
                productStuff[0]->lineData[i]));
        break;
    case UIntDS:
        for (int i = 0; i < width; i++)
            ((uint32_t*) fileData)[i] = round(
                productStuff[0]->calcOutputVal(
                productStuff[0]->lineData[i]));
        break;
    case FloatDS:
        for (int i = 0; i < width; i++)
            ((float*) fileData)[i] = productStuff[0]->lineData[i];
        break;
    case DoubleDS:
        for (int i = 0; i < width; i++)
            ((double*) fileData)[i] = productStuff[0]->lineData[i];
        break;
    default:
        printf("-E- OutFile_hdf4::writeLine - unrecognized data type = %d\n",
                productStuff[0]->dataStorage);
        exit(EXIT_FAILURE);
    }

    int32_t start[2];
    int32_t count[2];

    start[0] = currentLine;
    start[1] = 0;
    count[0] = 1;
    count[1] = width;

    if ((SDwritedata(sdsid, start, NULL, count, (VOIDP) fileData)) < 0) {
        printf("\n-E- OutFile_hdf4::writeLine(): SDwritedata unsuccessful\n");
        exit(EXIT_FAILURE);
    }

    if (qualityData) {
        if ((SDwritedata(quality_sdsid, start, NULL, count, (VOIDP) qualityData))
                < 0) {
            printf(
                    "\n-E- OutFile_hdf4::writeLine(): SDwritedata unsuccessful\n");
            exit(EXIT_FAILURE);
        }
    }

    currentLine++;
}

bool OutFile_hdf4::close() {

    if (fileData)
        free(fileData);

    if (metaData) {
        DPTB(
                SDsetattr(sdfid, "Data Bins", DFNT_INT32, 1, (VOIDP) & metaData->data_bins));
    }
    float tmpFloat = getFileMinVal();
    DPTB(SDsetattr(sdfid, "Data Minimum", DFNT_FLOAT32, 1, (VOIDP) & tmpFloat));
    tmpFloat = getFileMaxVal();
    DPTB(SDsetattr(sdfid, "Data Maximum", DFNT_FLOAT32, 1, (VOIDP) & tmpFloat));

    if (qualityData) {
        SDendaccess(quality_sdsid);
        quality_sdsid = -1;
    }
    SDendaccess(sdsid);
    sdsid = -1;
    SDend(sdfid);
    sdfid = -1;

    /*-----------------------------------------------------------------------*/
    /*  write map_palette */

    // make sure a palette has been loaded
    if (red && green && blue) {
        uint8_t data[768];
        int j = 0;
        for (int i = 0; i < 256; i++) {
            data[j++] = red[i];
            data[j++] = green[i];
            data[j++] = blue[i];
        }

        if ((DFPaddpal(fileName.c_str(), (VOIDP) data)) < 0) {
            printf("-E- OutFile_hdf4::close - Error writing map_palette.\n");
            return false;
        }

        uint16_t pal_ref;
        if ((pal_ref = DFPlastref()) > 0) {
            if ((DFANputlabel(fileName.c_str(), DFTAG_IP8, pal_ref,
                    (char*) "palette")) < 0) {
                printf(
                        "-E- OutFile_hdf4::close - Error writing palette label\n");
                return false;
            }
        }
    }
    return true;
}

int32_t OutFile_hdf4::addProduct(productInfo_t* productInfo) {
    return addProductNonDisplay(productInfo);
}

//------------------------------------------------------------------------------
// OutFile_netcdf4
//------------------------------------------------------------------------------

using namespace netCDF;

OutFile_netcdf4::OutFile_netcdf4() :
    OutFile() {
    ncFile = NULL;
    fileData = NULL;
    ncid = -1;
    qualityid = -1;
}

OutFile_netcdf4::~OutFile_netcdf4() {
    if (fileData)
        free(fileData);
}

void OutFile_netcdf4::setSize(int32_t width, int32_t height) {
    OutFile::setSize(width, height);
    if (fileData)
        free(fileData);
}

bool OutFile_netcdf4::open() {
    const char* tmpStr;
    currentLine = 0;

    try {

        // open file
        ncFile = new NcFile(fileName.c_str(), NcFile::replace);

        // write global metadata
        string prodName;
        size_t pos = fileName.find_last_of('/');
        if (pos == string::npos)
            prodName = fileName;
        else
            prodName = fileName.substr(pos + 1);
        ncFile->putAtt("product_name", prodName.c_str());

        ncFile->putAtt("instrument", metaData->sensor);
        ncFile->putAtt("title", metaData->title);
        ncFile->putAtt("project", PROJECT);

        if (strcmp(metaData->mission, "") != 0)
            ncFile->putAtt("platform", metaData->mission);

        ncFile->putAtt("temporal_range", metaData->prod_type);
        ncFile->putAtt("processing_version", metaData->pversion);

        strcpy(metaData->ptime, unix2isodate(time(NULL), 'G'));
        ncFile->putAtt("date_created", metaData->ptime);
        ncFile->putAtt("history", metaData->proc_con);
        ncFile->putAtt("l2_flag_names", metaData->flag_names);

        ncFile->putAtt("time_coverage_start",
                       unix2isodate(metaData->startTime, 'G'));
        ncFile->putAtt("time_coverage_end",
                       unix2isodate(metaData->endTime, 'G'));
        ncFile->putAtt("start_orbit_number", ncInt, metaData->start_orb);
        ncFile->putAtt("end_orbit_number", ncInt, metaData->end_orb);
        ncFile->putAtt("map_projection", mapProjection.c_str());

        ncFile->putAtt("latitude_units", metaData->lat_units);
        ncFile->putAtt("longitude_units", metaData->lon_units);

        ncFile->putAtt("northernmost_latitude", ncFloat, metaData->north);
        ncFile->putAtt("southernmost_latitude", ncFloat, metaData->south);
        ncFile->putAtt("westernmost_longitude", ncFloat, metaData->west);
        ncFile->putAtt("easternmost_longitude", ncFloat, metaData->east);
        ncFile->putAtt("geospatial_lat_max", ncFloat, metaData->north);
        ncFile->putAtt("geospatial_lat_min", ncFloat, metaData->south);
        ncFile->putAtt("geospatial_lon_max", ncFloat, metaData->east);
        ncFile->putAtt("geospatial_lon_min", ncFloat, metaData->west);
        ncFile->putAtt("grid_mapping_name", "latitude_longitude");

        float latStep = (metaData->north - metaData->south)
            / (float) getHeight();
        float lonStep;
        if (metaData->east < metaData->west)
            lonStep = (360 + metaData->east - metaData->west)
                / (float) getWidth();
        else
            lonStep = (metaData->east - metaData->west) / (float) getWidth();
        ncFile->putAtt("latitude_step", ncFloat, latStep);
        ncFile->putAtt("longitude_step", ncFloat, lonStep);
        ncFile->putAtt("sw_point_latitude", ncFloat,
                       metaData->south + latStep / 2.0);
        ncFile->putAtt("sw_point_longitude", ncFloat,
                       metaData->west + lonStep / 2.0);

        float geoResolution;
        const char* geoUnits;
        if (resolution > 1000.0) {
            geoResolution = resolution / 1000.0;
            geoUnits = "km";
        } else {
            geoResolution = resolution;
            geoUnits = "m";
        }
        char buf[2048];
        sprintf(buf, "%3.2f %s", geoResolution, geoUnits);
        ncFile->putAtt("geospatial_lon_resolution", ncFloat, geoResolution);
        ncFile->putAtt("geospatial_lat_resolution", ncFloat, geoResolution);
        ncFile->putAtt("geospatial_lat_units", metaData->lat_units);
        ncFile->putAtt("geospatial_lon_units", metaData->lon_units);
        ncFile->putAtt("spatialResolution", buf);

        ncFile->putAtt("number_of_lines", ncInt, height);
        ncFile->putAtt("number_of_columns", ncInt, width);
        ncFile->putAtt("measure", "Mean"); // FIXME

        ncFile->putAtt("suggested_image_scaling_minimum", ncFloat,
                       productStuff[0]->productInfo->displayMin);
        ncFile->putAtt("suggested_image_scaling_maximum", ncFloat,
                       productStuff[0]->productInfo->displayMax);
        if (!strcmp(productStuff[0]->productInfo->displayScale, "log"))
            tmpStr = "LOG";
        else if (!strcmp(productStuff[0]->productInfo->displayScale, "arctan"))
            tmpStr = "ATAN";
        else
            tmpStr = "LINEAR";
        ncFile->putAtt("suggested_image_scaling_type", tmpStr);
        ncFile->putAtt("suggested_image_scaling_applied", "No");

        ncFile->putAtt("_lastModified", metaData->ptime);
        ncFile->putAtt("Conventions", "CF-1.6 ACDD-1.3");
        ncFile->putAtt("institution", INSTITUTION);
        ncFile->putAtt("standard_name_vocabulary", STDNAME_VOCABULARY);
        ncFile->putAtt("naming_authority", NAMING_AUTHORITY);

        // create id
        strcpy(buf, metaData->pversion);
        if (strcmp(metaData->pversion, "Unspecified") != 0) {
            strcpy(buf, metaData->product_name);
            strcat(buf, "/L3/");
        } else {
            strcpy(buf, "L3/");
        }
        strcat(buf, metaData->product_name);
        ncFile->putAtt("id", buf);

        ncFile->putAtt("license", LICENSE);
        ncFile->putAtt("creator_name", CREATOR_NAME);
        ncFile->putAtt("publisher_name", PUBLISHER_NAME);
        ncFile->putAtt("creator_email", CREATOR_EMAIL);
        ncFile->putAtt("publisher_email", PUBLISHER_EMAIL);
        ncFile->putAtt("creator_url", CREATOR_URL);
        ncFile->putAtt("publisher_url", PUBLISHER_URL);

        ncFile->putAtt("processing_level", "L3 Mapped");
        ncFile->putAtt("cdm_data_type", "grid");

        // Some missions have DOIs
        if(clo_isSet(optionList, "suite")) {
            const char* doiStr = getDOI(metaData->mission, metaData->sensor, "L3M",
                                        clo_getString(optionList, "suite"),
                                        clo_getString(optionList, "pversion"));
            if (doiStr) {
                ncFile->putAtt("identifier_product_doi_authority", "http://dx.doi.org");
                ncFile->putAtt("identifier_product_doi", doiStr);
            }
            const char* keywordStr = getGCMDKeywords(clo_getString(optionList, "suite"));
            if(keywordStr) {
                ncFile->putAtt("keywords", keywordStr);
                ncFile->putAtt("keywords_vocabulary", KEYWORDS_VOCABULARY);
            }
        } // suite was set

        // create group: processing_control
        NcGroup grp1 = ncFile->addGroup("processing_control");
        grp1.putAtt("software_name", metaData->soft_name);
        grp1.putAtt("software_version", metaData->soft_ver);
        if ((tmpStr = strrchr(metaData->infiles, '/')) != NULL)
            tmpStr++;
        else
            tmpStr = metaData->infiles;
        grp1.putAtt("source", tmpStr);
        grp1.putAtt("l2_flag_names", metaData->flag_names);

        // create sub-group: input_parameters
        NcGroup grp2 = grp1.addGroup("input_parameters");
        char *end_str;
        char *token = strtok_r(metaData->input_parms, "|", &end_str);
        while (token != NULL) {
            char *end_token;
            strcpy(buf, token);
            char *name = strtok_r(token, "=", &end_token);
            for (uint32_t i = 0; i < strlen(name); i++) {
                if (name[i] == ' ') {
                    name[i] = 0;
                    break;
                }
            }
            strcpy(buf, strtok_r(NULL, "|", &end_token));
            grp2.putAtt(name, buf);
            token = strtok_r(NULL, "|", &end_str);
        }

        // define dimensions
        vector<NcDim> dimIds;
        dimIds.push_back(ncFile->addDim("lat", height));
        dimIds.push_back(ncFile->addDim("lon", width));

        /* (use old-style netCDF calls to define products) */
        ncid = ncFile->getId();
        int dimids[3];
        dimids[0] = dimIds[0].getId();
        dimids[1] = dimIds[1].getId();

        float tmpFloat;
        int status;
        idDS ds_id;
        ds_id.deflate = deflate;
        ds_id.fftype = DS_NCDF;
        ds_id.fid = ncid;

        size_t dataSize = 1;

        for (size_t prodNum = 0; prodNum < productStuff.size(); prodNum++) {

            int varDataType;
            switch (productStuff[prodNum]->dataStorage) {
            case ByteDS:
                varDataType = NC_BYTE;
                if (sizeof (signed char) > dataSize)
                    dataSize = sizeof (signed char);
                break;
            case UByteDS:
                varDataType = NC_UBYTE;
                if (sizeof (unsigned char) > dataSize)
                    dataSize = sizeof (unsigned char);
                break;
            case ShortDS:
                varDataType = NC_SHORT;
                if (sizeof (short) > dataSize)
                    dataSize = sizeof (short);
                break;
            case UShortDS:
                varDataType = NC_USHORT;
                if (sizeof (unsigned short) > dataSize)
                    dataSize = sizeof (unsigned short);
                break;
            case IntDS:
                varDataType = NC_INT;
                if (sizeof (int) > dataSize)
                    dataSize = sizeof (int);
                break;
            case UIntDS:
                varDataType = NC_UINT;
                if (sizeof (unsigned int) > dataSize)
                    dataSize = sizeof (unsigned int);
                break;
            case FloatDS:
                varDataType = NC_FLOAT;
                if (sizeof (float) > dataSize)
                    dataSize = sizeof (float);
                break;
            case DoubleDS:
                varDataType = NC_DOUBLE;
                if (sizeof (double) > dataSize)
                    dataSize = sizeof (double);
                break;
            default:
                cerr << "-E- OutFile_netcdf4::open - illegal data storage type = "
                     << productStuff[prodNum]->dataStorage << endl;
                exit(EXIT_FAILURE);
            }

            char* productNameStr = getProductNameFull(productStuff[prodNum]->productInfo);

            status = CreateNCDF(ds_id, productNameStr, /* short name */
                                productStuff[prodNum]->productInfo->description, /* long name */
                                productStuff[prodNum]->productInfo->standardName, /* NetCDF standard name (not set if passed NULL or "") */
                                productStuff[prodNum]->productInfo->reference,
                                productStuff[prodNum]->productInfo->comment,
                                productStuff[prodNum]->productInfo->units, /* units (not set if passed NULL or "") */
                                productStuff[prodNum]->productInfo->validMin, /* low end of valid range */
                                productStuff[prodNum]->productInfo->validMax, /* high end of range (no range set if low >= high) */
                                productStuff[prodNum]->productInfo->scaleFactor, /* scale factor (not set if 1.0)  */
                                productStuff[prodNum]->productInfo->addOffset, /* scaling offset (not set if 0.0)  */
                                productStuff[prodNum]->productInfo->fillValue, /* fill value */
                                varDataType, /* NCDF number type */
                                2, /* number of dimensions (must be <= 3) */
                                dimids /* dimension ids */
                                );
            check_err(status, __LINE__, __FILE__);

            int tmpVarId;
            status = nc_inq_varid(ncid, productNameStr, &tmpVarId);
            check_err(status, __LINE__, __FILE__);

            varid.push_back(tmpVarId);

            // add display info
            status = nc_put_att_text(ncid, tmpVarId, "display_scale",
                                     strlen(productStuff[prodNum]->productInfo->displayScale) + 1,
                                     productStuff[prodNum]->productInfo->displayScale);
            check_err(status, __LINE__, __FILE__);
            tmpFloat = productStuff[prodNum]->productInfo->displayMin;
            status = nc_put_att_float(ncid, tmpVarId, "display_min", NC_FLOAT, 1, &tmpFloat);
            check_err(status, __LINE__, __FILE__);
            tmpFloat = productStuff[prodNum]->productInfo->displayMax;
            status = nc_put_att_float(ncid, tmpVarId, "display_max", NC_FLOAT, 1, &tmpFloat);
            check_err(status, __LINE__, __FILE__);
        } // for prodNum

        // allocate the data for the biggest size
        if (fileData)
            free(fileData);
        fileData = allocateMemory(width * dataSize, "OutFile_netcdf4::open fileData");

        // add the quality variable
        if (qualityData) {
            productInfo_t* tmpInfo = allocateProductInfo();
            if (qualityName.empty()) {
                qualityName = (string) "qual_"
                    + getProductNameFull(productStuff[0]->productInfo);
            }
            status = findProductInfo(qualityName.c_str(), metaData->sensorID, tmpInfo);
            if (!status) {
                cerr << "-E- OutFile_netcdf4::open - cannot find product "
                     << qualityName.c_str() << " in product.xml" << endl;
                exit(EXIT_FAILURE);
            }

            status = CreateNCDF(ds_id, qualityName.c_str(), /* short name */
                                tmpInfo->description, /* long name */
                                tmpInfo->standardName, /* NetCDF standard name (not set if passed NULL or "") */
                                tmpInfo->reference, tmpInfo->comment, tmpInfo->units, /* units (not set if passed NULL or "") */
                                tmpInfo->validMin, /* low end of valid range */
                                tmpInfo->validMax, /* high end of range (no range set if low >= high) */
                                tmpInfo->scaleFactor, /* scale factor (not set if 1.0)  */
                                tmpInfo->addOffset, /* scaling offset (not set if 0.0)  */
                                tmpInfo->fillValue, /* fill value */
                                NC_BYTE, /* NCDF number type */
                                2, /* number of dimensions (must be <= 3) */
                                dimids /* dimension ids */
                                );
            check_err(status, __LINE__, __FILE__);

            status = nc_inq_varid(ncid, qualityName.c_str(), &qualityid);
            check_err(status, __LINE__, __FILE__);

            freeProductInfo(tmpInfo);
        }

        // Write lat/lon arrays if it is an SMI file
        if (mapProjection.compare("Equidistant Cylindrical") == 0) {
            float* lonarray = (float*) allocateMemory(getWidth()  * sizeof (float), "lonarray");
            float* latarray = (float*) allocateMemory(getHeight() * sizeof (float), "latarray");
            for (int i = 0; i < getWidth(); i++) {
                lonarray[i] = metaData->west + lonStep * i + lonStep / 2.0;
            }
            for (int i = 0; i < getHeight(); i++) {
                latarray[i] = metaData->north - latStep * i - latStep / 2.0;
            }

            productInfo_t *p_info;
            p_info = allocateProductInfo();
            NcVar var;

            // lat
            if (!findProductInfo("lat", metaData->sensorID, p_info)) {
                printf("lat not found in XML product table\n");
                exit(EXIT_FAILURE);
            }
            var = ncFile->addVar("lat", ncFloat, dimIds[0]);
            var.putAtt("long_name", p_info->description);
            var.putAtt("units", p_info->units);
            if (p_info->standardName != NULL) {
                var.putAtt("standard_name", p_info->standardName);
            }
            var.putAtt("_FillValue", ncFloat, p_info->fillValue);
            var.putAtt("valid_min", ncFloat, p_info->validMin);
            var.putAtt("valid_max", ncFloat, p_info->validMax);
            var.putVar(latarray);

            // lon
            if (!findProductInfo("lon", metaData->sensorID, p_info)) {
                printf("lon not found in XML product table\n");
                exit(EXIT_FAILURE);
            }
            var = ncFile->addVar("lon", ncFloat, dimIds[1]);
            var.putAtt("long_name", p_info->description);
            var.putAtt("units", p_info->units);
            if (p_info->standardName != NULL) {
                var.putAtt("standard_name", p_info->standardName);
            }
            var.putAtt("_FillValue", ncFloat, p_info->fillValue);
            var.putAtt("valid_min", ncFloat, p_info->validMin);
            var.putAtt("valid_max", ncFloat, p_info->validMax);
            var.putVar(lonarray);

            free(latarray);
            free(lonarray);
            freeProductInfo(p_info);
        }

        // store palette
        if (red && green && blue) {
            uint8_t data[768];
            int j = 0;
            for (int i = 0; i < 256; i++) {
                data[j++] = red[i];
                data[j++] = green[i];
                data[j++] = blue[i];
            }
            vector<NcDim> dimIds;
            dimIds.push_back(ncFile->addDim("rgb", 3));
            dimIds.push_back(ncFile->addDim("eightbitcolor", 256));
            NcVar var = ncFile->addVar("palette", ncUbyte, dimIds);
            var.putVar(data);
        }

    } catch (std::exception const & e) {
        cerr << "Exception: " << e.what() << endl;
    }

    return true;
}

void OutFile_netcdf4::writeLine() {
    int status;

    size_t start[3];
    size_t count[3];

    start[0] = currentLine;
    start[1] = 0;
    count[0] = 1;
    count[1] = width;

    for (size_t prodNum = 0; prodNum < productStuff.size(); prodNum++) {
        switch (productStuff[prodNum]->dataStorage) {
        case ByteDS:
            for (int i = 0; i < width; i++)
                ((int8_t*) fileData)[i] = round(
                    productStuff[prodNum]->calcOutputVal(
                    productStuff[prodNum]->lineData[i]));
            break;
        case UByteDS:
            for (int i = 0; i < width; i++)
                ((uint8_t*) fileData)[i] = round(
                    productStuff[prodNum]->calcOutputVal(
                    productStuff[prodNum]->lineData[i]));
            break;
        case ShortDS:
            for (int i = 0; i < width; i++)
                ((int16_t*) fileData)[i] = round(
                    productStuff[prodNum]->calcOutputVal(
                    productStuff[prodNum]->lineData[i]));
            break;
        case UShortDS:
            for (int i = 0; i < width; i++)
                ((uint16_t*) fileData)[i] = round(
                    productStuff[prodNum]->calcOutputVal(
                    productStuff[prodNum]->lineData[i]));
            break;
        case IntDS:
            for (int i = 0; i < width; i++)
                ((int32_t*) fileData)[i] = round(
                    productStuff[prodNum]->calcOutputVal(
                    productStuff[prodNum]->lineData[i]));
            break;
        case UIntDS:
            for (int i = 0; i < width; i++)
                ((uint32_t*) fileData)[i] = round(
                    productStuff[prodNum]->calcOutputVal(
                    productStuff[prodNum]->lineData[i]));
            break;
        case FloatDS:
            for (int i = 0; i < width; i++)
                ((float*) fileData)[i] = productStuff[prodNum]->lineData[i];
            break;
        case DoubleDS:
            for (int i = 0; i < width; i++)
                ((double*) fileData)[i] = productStuff[prodNum]->lineData[i];
            break;
        default:
            printf(
                    "-E- OutFile_hdf4::writeLine - unrecognized data type = %d\n",
                    productStuff[prodNum]->dataStorage);
            exit(EXIT_FAILURE);
        }

        status = nc_put_vara(ncid, varid[prodNum], start, count, fileData);
        check_err(status, __LINE__, __FILE__);
    } // for prodNum

    if (qualityData) {
        status = nc_put_vara(ncid, qualityid, start, count, qualityData);
        check_err(status, __LINE__, __FILE__);
    }

    currentLine++;
}

bool OutFile_netcdf4::close() {

    if (metaData) {
        ncFile->putAtt("data_bins", ncInt, metaData->data_bins);
    }
    ncFile->putAtt("data_minimum", ncFloat, getFileMinVal());
    ncFile->putAtt("data_maximum", ncFloat, getFileMaxVal());

    nc_close(ncid);
    ncid = -1;
    varid.clear();

    if (fileData) {
        free(fileData);
        fileData = NULL;
    }
    if (qualityData) {
        free(qualityData);
        qualityData = NULL;
    }

    return true;
}

int32_t OutFile_netcdf4::addProduct(productInfo_t* productInfo) {
    return addProductNonDisplay(productInfo);
}
